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1 INTRODUCTION 

Highly anisotropic layered superconductors can be considered as stacks of Josephson junctions (SJJ's) 
[0, 1^ . The key feature of SJJ's is mutual coupling of junctions, which can be achieved via magnetic 1^, Q , 
charging p| or non- equilibrium coupling mechanisms. Properties of SJJ's can be qualitatively 
different from that of single Josephson junctions (JJ's) or anisotropic type-II superconductors. This is 
reflected in a structure of magnetic vortices parallel and perpendicular to superconducting (S) layers. 
A perpendicular vortex represent an array of 2D pancake vortices in S-layers Those are Abricosov 
type vortices with normal cores and with supercurrents flowing in S-layers at a characteristic length Xat- 
Properties of such vortices are well studied, see e.g. Ref. ^, and will not be considered here. For 
parallel or slightly tilted magnetic field, vortex system commensurates with the layered structure [0 
due to intrinsic pinning which tends to locate vortices between S-layers. Such vortices (fluxons) are 
of Josephson type. They have circulating currents both along and across layers and do not have normal 
cores. 

Fluxon related properties of SJJ's were extensively studied both analytically and numerically [3,4,12- 
32], especially in connection with High-Tc superconductors (HTSC). Due to a high anisotropy, (weak 
Josephson coupling between layers), even small in-plane magnetic field causes penetration of fluxons. For 
Tl and Bi-based HTSC, the penetration field is ^ few Oe [33-35]. Since fluxon energy is low, they can be 
easily created by thermal excitations ^, |l^ or trapped during cooling down. Moreover, because fluxons 
can not move freely across layers, fluxon-antifluxon pairs in different junctions are stable and can exist 
even at zero magnetic field |3^ . Clearly, a knowledge of fluxon structure is necessary for understanding 
transport and magnetic properties of layered superconductors. From application point of view, it is 
particularly interesting to consider SJJ's with finite number layers, which is the case for HTSC mesas 
[30,36-43] and low-T^ SJJ's [30,44-54]. 

Here I will review anomalous single and multi- fluxon properties of SJJ's. Throughout the paper I 
will compare analytical, numerical and experimental results both for intrinsic HTSC and artificial low-Tc 
SJJ's. 

2 GENERAL RELATIONS 

The most common type of coupling in SJJ's is magnetic coupling. It appears when S-layers are thinner 
than London penetration depth, so that magnetic induction is not screened in one junction. SJJ's 
are coupled via shared magnetic induction and currents flowing along common electrodes. Properties 
of magnetically coupled SJJ's are described by a coupled sine-Gordon equation (CSGE), which was 
phenomenologically introduced in Ref. |^ and later rigorously derived in Ref. [^ 

First, I will briefly describe the formalism of CSGE, following notations of Ref. [Q. Let's consider 
a stack of N junctions with the following parameters: Jd -the critical current density, Ci and Ri- the 
capacitance and the quasiparticle resistance per unit area, respectively, ti- the thickness of tunnel barrier, 
di and Xsi - the thickness and London penetration depth of iS— layers, respectively, and L- the in-plane 
lengths of the stack. Elements of the stack are numerated from the bottom to the top, so that JJ i 
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consists of S-layers i, i + 1 and a tunnel barrier i. The CSGE can be written in a compact matrix form 



V?" = A • J, - Jb (1) 

Here <y9 is a column of gauge invariant phase differences ipi, 'primes' denote in-plane spatial derivatives, 
A is a symmetric tridiagonal N x N matrix with nonzero elements: Aij-i = —Si/Ai] Ai,i — Ai/Af, 

A-i^i+i = -Si+i/Ai, where + XsiCoth + \si+iCoth (^j^^ . '5'^ = Xs^cosech ■ Current 

density across layers, J^, consists of supercurrent, displacement current and normal current components: 

Jzi = jcisin {(fii) + CiCpi + OLi'^i. (2) 

Here jci = -j^, (7^ = 'dots' denote time derivatives, on = \J 2-Rc.hiCiii^ ^"^'^ damping parameters, 
$0 is the flux quantum and c is the velocity of light in vacuum. Jf, is a bias term [ po| . Space and 
time are normalized to Josephson penetration depth, Xji = sn^j'^iK ^^ ^^'^ inverse plasma frequency 

^/'^p' " respectively. 

Magnetic induction in a stack is given by: 

B = f A-V'. (3) 

Here = $o/(7rAj7A/). 

For d/Xs ^ 1, equations analogous to Eq.(l) were derived from Lawrence-Doniach (LD) ^ model 

fl The LD-CSGE conversion is shown in Table 1, along with estimations for Bi2212 HTSC (d=3 
, t =12 A, Jc=500-1000 A/cm2, A^=750-1000 A). 
In static case, CSGE reduces to 



ip" = A-jt,sin{(p). (4) 

For non-dissipative fluxon motion with a velocity u, the phase distribution remains unchanged in a 
coordinate frame moving along with the fluxon, ^ = a; — ut, and CSGE can be simplified ||32| 



D • 93" = A-jcSm(i^), (5) 



u^CA, and E is the unitary matrix. 
The first integral: CSGE has a first integral E 



-^'*Gp' +J2JC^[C0S{^^) - 1] = C, (6) 

where C is a constant of the first integral, (p'* is a string of (p[ and a matrix G = A^^ in static case and 
DA^^ for soliton-like fluxon motion. 

The free energy of the stack is a sum of magnetic, kinetic and Josephson energies. Using the first 
integral, Eq.(6), a simple expression for the free energy density is obtained in static case ||20|| : 



F = C + 2^j„[l-cos((^,)]. (7) 

Here F is measured in units of ^qJciXj /2Trc. From Eq.(7) it follows that the energy of any isolated 
solution (C ~ 0) is twice the Josephson energy, just like in single JJ. 



Table 1: Parameters of LD model in terms of CSGE and estimations for Bi2212 HTSC 
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Figure 1: Fluxon shapes in double SJJ's for different fluxon velocities: a) Solid and dotted lines represent 
numerical and analytic solutions, respectively, for ipi^2- b) Magnetic inductions Bi,2 obtained numerically. 
The existence of contracted and uncontracted components and the sign inversion of B2 (0) at u — > ci is 
clearly seen. From Ref. |25||. 



3 A SINGLE FLUXON 

Despite considerable theoretical efforts, there is still no exact fluxon solution for SJJ's. However, several 
approximate solutions were proposed. For infinite layered superconductor {N = 00) the first approximate 
solution was obtained within LD model ft was suggested that far from the fluxon center, current 
stream lines are elliptical with different London penetration lengths along and across the layers |fj] . The 
elliptic solution was extended by Clem and Coffey |l5| to explicitly take into account discreteness of 
layers and consider fluxon "core" region. Extension for superconductors with different layers was made 
in Ref. Et). The first approximate solution for finite SJJ's was derived for the case of two weakly coupled 
SJJ's [|l3t. However, since coupling strength was used as a perturbation parameter, that solution does 
not apply for strongly coupled SJJ's, which is the most interesting case. A "multi-component" solution, 
valid for arbitrary double SJJ's was obtained in Ref. pO| , using phase differences in fluxon-free junctions 
as a perturbation j25j. Such solution was shown to be in a good agreement with numerical simulations 
both for static and dynamic cases | |2^. Recently, the "multi-component" solution was extended for SJJ's 
with arbitrary number of junctions ]29|, |3^ . 

Apart from the multi-component, two other solutions were predicted in Ref. |^ for double SJJ's: 
(i) the single component solution, sin{ipi) ~ Ksin{ip2)', and (ii) a combination of the single component 
solution with a traveling plasma wave. Numerical modeling showed p6[ | that both solutions can be 
realized at high fluxon velocities. The latter solution was independently discovered in Ref. |Q and was 
interpreted as "Cherenkov radiation" from a rapidly moving fluxon. 

3.1 Approximate fluxon solution 

Let's consider arbitrary stack of N junctions with a single fluxon in junction zq. The problem with 
solving CSGE, Eq.(4), is coupling of nonlinear sin{ipi) terms in the right-hand side. Partial linearization 
with respect to ifi^ig was proposed in Refs.j2^, |32j for decoupling of CSGE. First Eq.(4) is diagonahzed 
to 
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Figure 2: Phase distribution for a fiuxon in a) N=ll and b) 600 SJJ's. Circles and lines represent 
numerical and analytic solutions, respectively. Insets show a) variation of magnetic induction and the 
derivative of phase in the central J J; b) variation of characteristic parameters of the fiuxon vs. the number 
of SJJ's. Data from Ref. || 



= sin{F^) + Er^, 



(8) 
(9) 



where = sin{ip,g) + J2 ''^m.Asin{ip,^,g) - sin{ip,„ + Y,i^m,iVi=tia)- Characteristic lengths, \^ and 

coefRcients ^ are given by eigenvalues and eigenvectors of A, respectively. The diagonalization proce- 
dure minimizes error functions, i?r„ far from the center. In this case Er^^ ~ ^1, while any other linear 
combination would yield ~ ipi. Therefore, Erm have a form of small ripple and vanish both inside and 
far from the fiuxon center. In the first approximation Er,-a terms in Eq.(8) can be neglected, leading to 
a set of decoupled sine-Gordon equations for Fm with well known solutions: 



F„i = iarctan (^e^/^'"^ . 
Finally, inverting Eq. (9) we obtain the approximate fiuxon solution p^, ISS 



(10) 



where K is the N x N matrix with elements equal to k 
given by 



ip = K-'F„,,iXrn), (11) 

For identical SJJ's eigenvalues/vectors are 



Am — Aj 



1 + 2Scos 



(-1) 



Tim 



-1/2 



iV+ 1^ 
i-io sin [-Kmi/ {N + 1)] 
sin [nmio/ {N -\- 1)] . 



m 1,2, . . .TV, 



(12) 
(13) 



Here S = Si/ Ki is the coupling parameter, S" = 0.5. Normalization of eigenvectors in Eq.(13) follows 
from Eq.(9) and ensures that the total phase shift is 27r for i = io and zero for the rest of the junctions. 

For SJJ's with odd and a fiuxon in the middle junction, — A^/2, the number of components reduces 
to n = {N + l)/2, due to a symmetry relation, (pi„-j = y^io+j^ ^^'^ solution becomes particularly 
simple: 



4 




Figure 3: Numerically simulated distributions of magnetic induction for the fluxon in iV=600 SJJ's: a) 
B{z — const) along the a6-plane. Dashed line shows the derivative of phase in the central junction. 
Inset represents comparison between numerical, multi-component and elliptic solutions in the central 
junction, b) B{x = const) along the c-axis. Inset represents a 3D plot of B{x, z) in the fluxon core 
region. From Ref. 



f^o = (TO = l,3,...iV), (14) 



Here are examples of fluxon solutions for small N: 



N ^2 



<P2 -~ 



Ei±E^,\^^\j{l + S)-^''' 
\j{l^S)~^'^ 



2 '^^2 



^1,3 = A3 -A, (1-725)-'^ 



f ^3-™^,Al=A,;(l + ^/35)-'/' 

A characteristic feature of the multi-component solution is that the fluxon is characterized by N 
lengths, Am, different from the two lengths of the in elliptic solution ||l^, |l^, |l5|. The approximate fluxon 
solution for the dynamic CSGE, Eq.(5), is given by the same expression, Eq.(ll), but with renormalized 
characteristic lengths. 



Cm, 



(16) 



with characteristic velocities 
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Figure 4: Numerically simulated current-velocity characteristics for a single JJ iV = 1 (dashed line), and 
for SJJ's with different TV: TV = 2 (sohd hue), iV = 3 (open circles), iV = 5 (dotted hne) and iV = 21 
(filled triangles). From Ref-H 
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, m = 1,2...7V. (17) 



Here cq = c-^/i/e^A is Swihart velocity for a single junction. 

In Fig. 1 spatial distributions of a) phase and b) magnetic induction are shown for a fluxon in two 
strongly coupled (5 = 0.495) identical SJJ's at different velocities. Parameters of the stack are: di-s = 
t\,2 = O.OlAj, Asi-3 = O.lAj. A quantitative agreement between "exact" and approximate solutions is 
seen. For identical double SJJ's exactly one half of the fluxon belongs to each of the components, F\,2- 
Indeed, from Fig. 1 a) it is seen that for w ~ ci there is a Lorentz contracted core at x = with one-7r 
step in On both sides of the core, there are two 7r/2 tails, decaying at distances ^ \\. Fig. 1 

b) represent numerically simulated profiles of i?i^2(a;). It is seen that the fluxon shape becomes unusual 
in dynamic state: a dip in Bi develops in the center of the fluxon with increasing velocity and at m — > ci 
-62(0) changes sign. At u = ci, i?2(0) — — i?i(0), in agreement with the multi-component solution. 

Numerically simulated fluxon profiles for N —W and 600 SJJ's are shown in Figs. 2 and 3. Parameters 
are typical for Bi2212: d==3 A, t =12 A, Aj = l/im (J^ ~ GTSA/cm^), As=750 A(Aafc ~ 0.17^to). It was 
shown that the multi-component solution correctly describes the fluxon shape in static case for arbitrary 

sj J 01111,11. 

3.2 Decoupling of phase and field 

Inset in Fig. 2 a) shows spatial variation of B and Lp' in the central JJ zo=6 of TV = 11 SJJ's. It 
is seen that length scales for variation of (/? and B are different: varies at distances ~ Aj, while the 
characteristic length for variation of B at large distances, A(oo), is An — 5.4Aj, Eq.(12). The discrepancy 
between Lp\^ and becomes dramatic for large A^, see Fig. 3 a). This is in sharp contrast with behavior 
of Josephson vortices in single JJ's or Abricosov vortices in type-II superconductors, in which spatial 
variations of both phase/current and magnetic induction are given by the same length scale. Existence 
of one and the same length scale for / and B is viewed as a natural consequence of Ampere's law. 

So, how could phase/current decouple from magnetic induction in SJJ's? This can be understood 
from the multi-component fluxon solution: (i) Magnetic induction in SJJ's is non-local, Bi depends on 
kpi in all junctions, see Eq.(3), which is an essence of magnetic coupling, (ii) If the fluxon has a multi- 
component structure, then distribution of phase is determined by the shortest, while magnetic induction 
- by the longest characteristic length. Indeed, according to Eq. (14), the effective Josephson penetration 
depth in the center of the fluxon is 
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Figure 5: Phase distributions for the fluxon moving with velocity u ~ 0.928co in = 11 SJJ's are shown 
for the central junction, i = 6, (top panel); i = A (middle panel); and i = \ (bottom panel). Spectra of 
Cherenkov radiation are shown in insets a) for i = 6 and b) for « = 4 (solid line) and i = 1 (dotted line) . 
FromRef. ||. 
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which is obviously dominated by the shortest length Ai ~ Aj. On the other hand, from Eq. (15), variation 
of B is dominated by the longest length, Ajv, which according to Eq. (12) increases rapidly with TV for 
strongly coupled SJJ's, S ^ 0.5. 



3.3 N-dependence and assymptotics 

Fluxon shape strongly depends on the number of junctions in a stack. Inset in Fig. 2 b) summarizes 
variation of the fluxon shape for different N . Open circles represent magnetic induction at the center of 
the fluxon -B(O). -B(O) increases with N and saturates at i?(0) ~ 4iJo for N > Xab/s. Triangles represent 
numerically obtained A(0), which characterize variation of phase/current in the fluxon core. In agreement 
with Eq. (18), A(0) increases only slightly with and saturates at ~ 1.3A,/ for > 10. On the other 
hand, the effective magnetic length far from the core, A(oo), increases dramatically with N. A(cxd) is 
determined by the largest characteristic length Aat, as shown by the solid line. For large N > Xab/s, 
Xn ^ Ac, see dotted line in inset of Fig. 2 b), and the multi-component solution approaches the elliptic 
solutions |l^. This is demonstrated in inset of Fig. 3 a). Top axes in Figs. 3 a) and b) show that 
B{x, z) varies at length scales ~ Ac and A^h along and across layers, respectively. It is also seen that the 
most spectacular feature of the fluxon is a sharp peak B{Q,0), representing the "Josephson core" of the 
fluxon '--^ Aj X s in a6-plane and c-axis, respectively, as shown in inset of Fig. 3 b). 

Since the fluxon has several characteristic lengths, a question may arise, what is the effective Josephson 
penetration depth, Aj(e^j). This is important to know, because behavior of SJJ's changes drastically 
when L > Xj[f.ff)- Roughly speaking, fluxons can exist only in "long" JJ's, L 3> -^./(e//): causinga large 
difference in behavior of "short" and "long" JJ's. Length scales somewhat different from that in JlJ] were 
derived in Ref. ||lj], based on simulations for N — 15. In general, different conclusions about the fluxon 
size could be made from numerical simulations for different N: short scale ^ Xj for N=7 [Q, and A^=19 
[ p^ ; and long scale ^ Ac for A^=50 Long scale, ^ Ac, variation of B for N > 1000, was deduced 

from experimental observations of the fluxon in HTSC [Q. From the analysis above it is clear that such 
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Figure 6: Gray scale plot of Gibbs free energy versus the number of fluxons in double SJJ's a± H = l.lHci- 
Darker regions correspond to smaller Gibbs energy. The existence of multiple quasi-equilibrium states is 
seen. From Ref. IM. 



"discrepancy" is a consequence of strong A''— dependence of the fluxon shape, however, Aj(eyy) = A(0) is 
always — Aj, as shown in inset of Fig. 2 b). 



3.4 Flux-flow characteristics 

Variation of fluxon shape in dynamic state is particularly interesting both because it has a direct impact 
on c— axis transport properties and because fluxon shape becomes very unusual at high propagation 
velocities, see Fig. 1. When bias current is applied, a fluxon starts moving along the junction with a 
velocity determined by the balance between Lorentz and viscous friction forces. This results in appearance 
of a flux-flow branch in current-voltage characteristics (JVC's) with average DC voltage Vig/Vo = u/cq, 
{Vq — -^If^) in the JJ containing the fluxon and zero for i ^ iq. Therefore, current-velocity characteristics 
are equivalent to IVC's. In a single JJ, fluxon is a relativistic object. It gets Lorentz contracted |Q when 
the fluxon velocity approaches Swihart velocity, i.e. the phase velocity of electromagnetic waves in the 
transmission line formed by the junction. However, in SJJ's according to the multi-component solution, 
only Fi component of the fluxon gets Lorentz contracted at the lowest velocity ci . The weight of this 
component decreases with iV, i.e., there is less contraction for larger N . 

So, what happens with the fluxon in SJJ's: does it gets contracted at u = ci or not? Systematic 
analysis of the double SJJ's showed that the answer strongly depends on parameters of SJJ's Q. From 
Fig. 1 it is seen that the fluxon in identical double SJJ's is well described by the multi-component 
solution up to u = Ci. However, in general, partial linearization procedure fail at u = Ci because (pi^i^^ 
no longer small. It was found [^, that for > 1 and K2 > 1, a "single component" (F2) solution [go) 

is realized in double SJJ's. This solution was anticipated in Ref. and is characterized by relationship 
sin{Lpi) ~ K2sin{ip2)- Note, that it becomes exact at u = ci for the above mentioned parameters. 
Apparently, the single F2 component fluxon does not contract at m = ci. 

In Fig. 4, numerically simulated IVC's are shown for a single junction = 1, and SJJ's with TV —2, 
3, 5 and 21. Simulations were done for annular junction geometry with L = lOOAj, parameters typical 
for Bi2212 HTSC and damping coefficient ai — 0.05. Vertical dotted lines represent lowest characteristic 
velocities: ci{N = 2) ~ 0.8165, ci{N = 3) ~ 0.7653 and ci{N = 21) ~ 0.7071co. Clear velocity matching 
step at u ^ cq is seen for a single junction Vertical portions of IVC's at u ^ Ci are also seen for 
N — 2 and 3, as a result of partial Lorentz contraction of the fluxon, see Fig. 1. However, for iV = 5 
and 21 there is no peculiarity in IVC's at a velocity matching condition u = ci, indicating absence of 
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Figure 7: Three simulated /c(-ff) sub-branches of double SJJ's, corresponding to submodes of the (2,2) 
mode with different fluxon sequence. Insets a-c) show spatial current distributions in junctions 1 and 2 
for each submode at i?/iJo = 80. From Ref. |2§ 



fluxon contraction. In Rcf. [ p2| it was found out that Lorcntz contraction at u = ci is absent for SJJ's 
with iV > 3. Note that there are very little changes in IVC's for > 5 (the IVC's for TV = 5 and 21 
are almost indistinguishable in Fig. 4). This is due to the fact that dissipation is mostly determined by 
the fluxon core region, which remains almost unchanged for N > 5^ see A^— dependence of A(0) in inset 
of Fig. 2 b). 

3.5 Cherenkov radiation 

Due of the lack of Lorentz contraction, the fluxon can propagate with a velocity larger than ci as 
independently proposed in Refs. ||5^, Such super-relativistic motion is accompanied by "Cherenkov" 
radiation jl^, ^ ^ . An important insight to Cherenkov radiation can be obtained from the 

multi-component solution. For u > c„, characteristic lengths Xm<n become imaginary, characteristic for 
"wave" solutions. Indeed, far from the fluxon center we can further linearize Eq.(8), and obtain that the 
solution for Fm,<n is a traveling wave with the dispersion relation ||20|, |3^: 



l^7n = "^/'^p - ^mJ ^ = UKm, (19) 

where uj^ = Wpy^l — J^- Those are Josephson plasma waves Q moving along with the fluxon at 
velocity u p^ . 

In Fig. 5, numerically simulated phase distributions for a rapidly moving fluxon (cg > u ~ 0.928co > 
C5) in iV = 11 SJJ's are shown for a central junction, containing the fluxon, — Q (top panel); i — A 
(middle panel); and the outmost junction i = 1 (bottom panel). Parameters of SJJ's are the same as in 
Fig. 4. The six lowest characteristic velocities for iV = 11 are ci = 0.7132co, C2 — 0.7320co, C3 — 0.7653co, 
C4 — 0.8165co, C5 = 0.8913co, cg = cq. Well defined oscillations are seen behind the fluxon (the fluxon is 
moving from right to left). The spectra of oscillations are shown in insets a) for z = 6 and b) for i — A 
(solid line) and i = 1 (dotted line). Clear beatings of oscillations are seen in JJ's 6 and 1, which is a 
signature of coexistence of multiple wave lengths. Indeed, spectra of oscillations exhibit several maxima 
corresponding to odd plasma modes, as indicated by arrows in insets a and b). 

Numerical simulations support the assumption po| that Cherenkov radiation at u > c„ is due to 
degeneration of fluxon components Fm<n into Josephson plasma waves ||2^, Namely, it was observed 
that: (i) Plasma mode m is generated when u > Cm- (ii) Only those modes which constitute the fluxon, 
appear in Cherenkov radiation. E.g., for a fluxon in the central junction of a stack with odd N there 
are no even modes, see Eq. (14) and Fig. 5. (iii) Relative amplitudes of oscillations are related to 
weights coefficients of the corresponding fluxon component. For the case of Fig. 5, weight coefficient of 
mode 3 in junction 4, = 0, and plasma mode 3 is absent in the spectrum. On the contrary, for the 
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Figure 8: a) A gray scale plot of the probability distribution of Ic{T) for a Bi2212 mesa. Inset demon- 
strates multiple-peak structure of P (Ic) for several T. b) Simulated Ic{T) for three SJJ's with Bi2212 
parameters at H=0. Multiple Ic{T) branches are clearly seen. Branches corresponding to (0,0,0) and 
simplest fluxon-antifluxon modes are marked. Insets show distribution of fluxons in the stack. From Ref. 
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outmost junction, the weight coefficient of mode 3, 3 ~ 0.118, is significantly larger than for mode 1, 
K^l ~ —0.043, and the maximum of amplitude corresponds to plasma mode 3, see inset b) in Fig. 5. 

4 MULTI-FLUXON CONFIGURATIONS 

Due to repulsive fluxon interaction, fluxons arc believed to form a regular fluxon lattice (FL) at large 
parallel magnetic field. Since magnetic induction of a fluxon is spread over many layers, see Fig. 3, it 
is not favorable to have a perfect translational symmetry along the c— axis, and FL is expected to be 
rhombic [ ]60[ . However, the perfect FL can be achieved only when fluxons strongly interact with each 
other, i.e., when spacing between fluxons along the a6-plane is ~ Aj. For Bi and Tl based HTSC, the 
corresponding magnetic field, <I>o/A,/s ~ IT is three orders of magnitude larger than the lower critical 
field iyl'i -mT H, g. 

Fluxon distribution in layered superconductors at low fields, h]}^ < H < ^q/Xjs, is still a matter 
of controversy. Levitov |Q have found " phyllotaxis" and bifurcations in the FL at low fields, due to 
existence of multiple metastable FL's with approximately equal energies. In Ref. ]6l| |, this approach 
was extended using a "growth" algorithm, allowing variation of the FL along the c-axis. Although 
restricted to constant space periodicity along layers, the model showed the existence of a large variety of 
quasi-periodic or aperiodic fluxon configurations in the c-axis direction. At low fields, the fluxon system 
becomes completely frustrated and tends to be chaotic . Recent numerical simulations revealed that 
fluxons in SJJ's may form buckled chains instead of the regular FL. Clearly, the lattice description 

of multi- fluxon distribution becomes unappropriate at low magnetic fields. 

4.1 Quasi-equilibrium fluxon modes and submodes 

Recently it was shown that multiple quasi-equilibrium fluxon configurations (modes) exist in long, 
L '3> \j, strongly coupled SJJ's |Q. Fig. 6 shows a numerically simulated gray scale plot of Gibbs free 
energy versus the number of fluxons in double SJJ's for applied parallel field slighltly above the lower 
critical field, H — l.lHd- Parameters of the stack are Jc2 = 2Jci, ^Si = O.lAji, ti = di = O.OlAji, 
L = 50Aji. Darker regions correspond to smaller Gibbs energy. The absolute minimum, representing 
a thermodynamic equilibrium, correspond to 7 fluxons in junction 1 and zero fluxons in junction 2. 
However, it is seen that there is a number of other fluxon configurations (modes) for which local minima 
of Gibbs free energy is achieved. All those modes are stable and represent quasi-equilibrium states of 
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Figure 9: a) /c (T), for Nb/Cu multilayer. Different 
strong fluctuations of Ic take place in the 2D state, 
and perpendicular (open circles) to layers. A kink 
Ref. [|0| 



symbols correspond to different runs. It is seen that 
b) The upper critical field parallel (solid rectangles) 
in represents the 3D-2D crossover. Data from 



the stack. The modes are characterized by different fluxon distributions in the stack and do not have 
a —priori any translational symmetry. We will use a string (ni, ...ni, ...n^) as a short notation of fluxon 
modes, where rii is a number of fluxons in junction i and N is the number of junctions in the stack. E.g., 
the equilibrium state in Fig. 6 corresponds to (7,0) mode. 

In Ref. p3| existence of fluxon submodes was demonstrated, which have the same number of fluxons 
but are different with respect to fluxon sequences and symmetry of phase differences in the junctions 
(along the planes). Insets a-c) in Fig. 7, show three different submodes of (2,2) mode for a double stack 
with the same parameters as in Fig. 6. Those submodes have nearly equal self energies, i.e. nearly equal 
probability, however, correspond to distinctly different Ic{H) sub-branches. 

The number of fluxon modes and submodes increases rapidly with increasing number of junctions and 
fluxons. From mathematical statistics it follows that for a stack with N non-identical junctions with M 
fluxons, the total number of fluxon modes (indistinguishable fluxon sequences) is 



_ (N + M-iy. 

modes - , l^Uj 



and the number of submodes (distinguishable fluxon sequences) is 



^submodes ■ (^-^) 

A specific feature of SJJ's is existence of stable fiuxon-antifluxon pairs |Q. Due to attractive in- 
teraction between fluxons and antifiuxons, they tend to collapse both in a single junction or type-II 
superconductor. However, in SJJ's fiuxons can not move across S-layers. Therefore, a fluxon and an 
antifluxon in neighbor junctions do not collapse, but form a bound pair. Furthermore, once in the stack, 
such pair can not be removed (unless destroyed) by transport current because there is no net Lorentz 
force on the pair. Fiuxon-antifluxon pairs are likely to be present in SJJ's at H < Hd and can dominate 
transport properties at zero field |3^, |3l) . 

Due to existence of multiple quasi-equilibrium fiuxon (and fiuxon-antifluxon) modes/ submodes, the 
state of SJJ's is not well defined. It can be described only statistically with a certain probability of being 
in any of the modes. Neither the number of fiuxons in each junction nor even the total number of fiuxons 
in the stack is fixed for given H and T. Moreover, since different fluxon modes /submodes may have 
similar energies, the system is frustrated and exhibit strong fluctuations. This may account for unusual 
c-axis transport properties of both HTSC mesas ^ and low-Tc SJJ's |5^, |30t| . 
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Figure 10: Ic{H) dependencies for a) Nb/Cu multilayer, and b) Bi2212 mesa. Symbols correspond to 
maxima in P{Ic), shown in inset. Data from Ref. |]30| 



4.2 Multiple-valued critical current 

Statistical analysis is crucial for probing multiple quasi-equilibrium states in SJJ's. Fig. 8 a) presents 
in a gray scale a probability distribution of critical current P (/c) vs. T for 50 fim long Bi2212 mesa 
with N — 5 intrinsic SJJ's, at earth magnetic field. Darker regions correspond to a larg er p robability. 
To obtain P (Ic), 5120 switching events from zero- voltage state were measured at each T j62[ Inset 
in Fig. 8 a) shows P {Ic) for 40 K < T < 60 K. It is seen that P (Ic) is very wide and consists of several 
superimposed peaks, attributed to different fluxon modes 

In Fig. 8 b) simulated Ic{T) dependencies are shown for = 3 long (L ~ lOAj) SJJ's with Bi2212 
parameters, at H~0. Different symbols represent separate runs with increasing or decreasing T and/or 
different initial conditions. It is seen that Ic{T) consists of multiple distinct branches, in qualitative 
agreement with experiment. Insets in Fig. 8 b) show amount of fluxons in the junctions. It can be 
seen that each Ic{T) branch corresponds to a particular fluxon or fluxon- antifluxon mode. I want to 
emphasize that Ic is multiple valued even at H — 0. This is due to stability of fluxon-antifluxon pairs in 
SJJ's as discussed above. Branches corresponding to simplest fluxon-antifluxon modes (-1,1,0), (0,-1,1) 
and (1,0,-1) are marked in Fig. 8 b). 

Fig. 9 a) shows Ic (T) for a Nb/Cu (20/15 nm — 10) multilayer. Different symbols correspond to 
different runs. Measurements were done at ^ 10 mT along layers. It is seen that Ic strongly fluctuates 
below ~5 K. Fig. 9 b) shows the upper critical field IIc2{T) parallel (solid rectangles) and perpendicular 
(open circles) to layers. A kink in 7f|2(r) reflects dimensional 3D-2D crossover, which is well studied 
in those systems [48,53,63-66]. The 3D-2D crossover is one of pecuhar properties of superconducting 
multilayers ||67|. In the 3D state a multilayer behaves as a bulk superconductor, while in the 2D state 



- as a stack of Josephson junctions 53 . Dashed lines in Fig. 9 show, that strong fluctuations of Ic 
appear in the 2D state reflecting appearance of SJJ's. This unambiguously demonstrates that multiple 
valued critical current is a genuine feature of SJJ's, due to existence of multiple fluxon modes. 

4.3 Magnetic field dependence of the critical current 

The dependence of Ic on in-plane magnetic fleld is a crucial test for dc Josephson effect. In a single J J, 
Fraunhofer modulation of IdH) occurs as a resvflt of flux quantization [^. For SJJ's it was predicted 
that the periodicity of IdH) is ||6^ 



AH = ^o/Ls. 



However, experimental Ic{H) both for HTSC ||70|, |7lJ g3[ and low-Tc 
such periodicity. So, what's wrong with dc- Josephson effect in SJJ's? 



(22) 

19[ 0, H SJJ's do not show 
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Figure 11: Simulated Ic{H) for a) long, strongly coupled b) long, weakly coupled and c) short, strongly 
coupled double SJJ's. It is seen that Fraunhofer modulation is obscured in long, strongly coupled SJJ's 
due to existence of multiple quasiequilibrium fluxon modes. However, periodicity of Ic{H) is restored 
either in weakly coupled or short SJJ's. Data from Ref. |E3| 



Experimental Ic{H) for Nb/Cu multilayer (A^=10, L =20 ^ m) and Bi2212 mesa {N =5, L = 20 ^m) 
are shown in Fig. 10 a) and b), respectively, at T = 4.2 K. Note a logarithmic scale of Ic- In Fig. 10 a), 
different symbols represent different zero- field-cooled runs with increasing field. Fig. 10 b) shows a result 
of statistical analysis. Large and small symbols in Fig. 10 b) represent main and secondary peaks in 
P(/c), respectively, for increasing (circles) and decreasing (diamonds) field. The probability distribution 
with several maxima is shown in inset. It is seen that Ic is multiple valued and, even though one could 
probably recognize certain branches, there is no periodic Fraunhofer modulation. 

Fig. 11 presents numerically simulated Ic{H) patterns for a) long strongly coupled {L = lOAji, 
S = 0.495, d/Xs = 0.1); b) long weakly coupled (L = lOAji, S = 0.0064, d/Xg = 5) and c) short strongly 
coupled {L = 2Xji, S = 0.495, d/Xs — 0.1) double SJJ's with Jc2 — 2Jci- Note logarithmic scale for Ic- 
Plots contain several runs for increasing or decreasing field with different initial conditions. From Fig. 
11 a) it is seen that the Ic{H) pattern for long, strongly coupled SJJ's consists of multiple closely spaced 
branches. Each branch is characterized by a certain fluxon distribution, i.e., corresponds to a certain 
fluxon mode or submodes, see Fig. 7. The Ic{H) pattern does not exhibit clear periodicity because Ic 
switches randomly between closely spaced branches. This does not mean that something is wrong with 
dc Josephson effect, but rather is a consequence of magnetic coupling between junctions. 

From Figs. 11 b) and c) it is seen that periodicity of IdH) is restored either when the coupling 
is decreased or a stack becomes short compared to A,/. For weakly coupled SJJ's, Ic{H) is simply 
determined by the smallest critical current of individual JJ's, while in short SJJ's there are no fluxons 
and, consequently, no fluxon modes. A transition from aperiodic to periodic modulation of IdH) with 
decreasing L/Xj was observed both for Nb/Cu multilayers |5^ and for YBa2Cu307_a; |^ close to Tc- 

4.4 Multiple flux-flow sub-branches 

In sec. 3.4 we considered viscous motion of a single fluxon. For a multi-fluxon configuration, time 
averaged flux- flow voltage is proportional to the number of fluxons in SJJ's and a fluxon velocity. For 
large H, the ratio V/H is: 



V/H ~ Nsu. (23) 

Note that it is independent on the junction length. A characteristic feature of SJJ's is splitting of Swihart 
velocities, see Eq.(17), and corresponding splitting of flux-flow IVC's [|[ ^ Flux-flow phenomenon 
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Figure 12: a) Flux-flow IVC's of 50 long Bi2212 mesa with TV = 5 intrinsic SJJ's at r=4.2 K are 
shown for H from 0.148 to 0.396 T with a step ~0.025 T. b) Enlarged top parts of / vs. V/I curves 
are shown for H varying with a small step ~5 mT in the range iJ=0. 351-^0. 396 T, at T=4.2 K. Closely 
spaced flux-flow sub-branches are seen. Data from Ref. ||30] 



was intensively studied both for HTSC [30,37-41,74,75] and low-T^ © 1^, H SJJ's due to possible 
application as a source of microwave radiation 1 27 . From application point of view, the most interesting 
is coherent (in-phase) fluxon propagation Q at the highest characteristic velocity cat, because it could 
generate narrow linewidth radiation with large amplitude |2^. However, so far the measured radiation 
linewidth from SJJ's was much wider than expected ||7^, [55|. 

In Fig. 12 a) flux-flow IVC's of long [L = SO^im, N = b) Bi2212 mesa are shown for H = 0.148-0.396 
T with a step ~ 0.025 T, at T=4.2 K. Magnetic fleld was applied along the a6-plane and perpendicular 
to the longest side of the mesa. It is seen, that a low resistance branch develops in IVC's with increasing 
H . The 20 /im long mesa on the same single crystal showed the same V/ H value, in agreement with Eq. 
(23). However, IVC's were more complicated due to presence of Fiske steps The maximum flux-flow 
velocity is close to the lowest characteristic velocity ci ~ 2.5 — 3 x \{fcm/s. This value is consistent with 
observations by different groups [37-41,74,75] and is in agreement with the Fiske step voltage observed 
in Bi2212 mesas [0. 

From Fig. 12 a) it is seen that flux-flow branches exhibit strong fluctuations. At small voltages, the 
spread in IVC's reflects fluctuations in Ic, as discussed in sec. 4.2. In Fig. 12 b) enlarged top parts of / 
vs. R = V/I curves for the same mesa are shown at magnetic fields in the range 0.351 - 0.396 T, varying 
in small steps ~ 5 mT, at T=4.2 K. The characteristic feature of IVC's from Fig. 12 is that the curves 
are not just uniformly wide, but rather consist of multiple closely spaced but distinct sub-branches. The 
IVC's switch hysteretically between sub-branches when current is swept back and forth. Varying i?, we 
only change the set of available sub-branches but do not change the shape of a particular sub-branch. 

Similar multiple flux-flow sub- branches were observed for Nb/Cu multilayers [Q and HTSC mesas 
[[74[ [ttI . In the latter two papers it was suggested that sub-branching is due to splitting of Swihart 
velocities. Indeed such splitting was observed for low-T^ SJJ's Q and possibly for HTSC mesas 

[ [75| . In those cases branches revealed different V/H scaling, corresponding to different c„, see Eq.(23). 
On the contrary, tiny sub-branches shown in Fig. 12 b) (as well as those observed in Refs. [^ have 
approximately the same V/H scaling corresponding to the lowest velocity ci. Furthermore, because of 
small N=5, splitting of c„ should correspond to almost two orders of magnitude larger splitting than 
that in Fig. 12 b) [^. It should also be noted that the number of sub-branches is not limited by the 
number of SJJ's in the mesa. All this rules out explanation of multiple flux- flow sub-branches, shown in 
Fig. 12 b), in terms of splitting of Swihart velocities. 

In Refs. [m| it was proposed that multiple flux-flow sub-branches correspond to propagation of 



different quasi-equilibrium fluxon modes. This was supported by direct numerical simulations [Eg, 30 



Fig. 13 presents numerical modeling of IVC's for the mesa from Fig. 12 ( iV = 5, Jc ~ 10^/cm 
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Figure 13: Simulated IVC's are shown for two back-and-forth current sweeps (denoted as A and B) for 
the Bi2212 mesa from Fig. 12. Top three panels show overall IVC's V — J^^i^ with successive "zoom- 
in" to a small voltage region from panel a) to c). Fluxon modes realized at the end and the beginning 
of flux-flow branches for both sweeps are specified. In the bottom row, IVC's of individual junctions 
(successively offset along the V— axis for clarity) and the number of fluxons in the junctions are shown. 
It is seen that fluctuation of the flux-flow voltage is due to switching between fluxon modes. 



s = 15 A, d — 3A, A J ~ 1/im, L/Xj = 50 and damping parameters ai — 0.05). Junctions are identical, 
except for larger thickness of the bottom electrode, di — 6A, in order to imitate the bulk crystal beneath 
the mesa. Two successive back-and-forth current sweeps are shown (denoted as A and B) at in-plane 
magnetic field H = W^Hq ~ 2000 Oe. Top three panels show overall IVC's of the stack, V = J^^i^ 
with successive "zoom- in" to a small voltage region: panel a) demonstrates IVC's in a large scale, panel 
b) presents the flux-flow branch, and panel c) shows details of switching between supercurrent {V = 0) 
and flux-flow branches. In the bottom row, the left and right groups of plots represent data for sweeps 
A and B, respectively. In each group, the left plots represent individual IVC's of all five junctions 
(successively shifted along V— axis for clarity), while the right plots show amount of fluxons in each 
junction, determined as [ipi{L) — ipi{0)]/2T: (lines), and the total amount of fluxons in the stack (symbols). 
Despite nominally the same external conditions, the IVC's A and B are apparently different. For the 
IVC-A only three (z =2, 3, and 4), while for the IVC-B all five junctions switch to the normal branch 
(shown by dotted lines in Fig. a)) at the end of the flux-flow branch. 

From Fig. 13 b) it is seen that flux-flow branches exhibit strong fluctuations, caused by variation of 
fluxon numbers in junctions, as shown in the bottom row of Fig. 13. Each time fluxons leave the stack on 
one side, other fluxons entering the stack on the opposite side have freedom to choose between all possible 
quasi-equilibrium fluxon modes. Each fluxon mode is represented by a distinct flux-flow sub-branch, with 
slightly different Ic and flux- flow voltage. Two of such sub- branches can be seen at I/Ic = 0.1 — 0.15 in 
Fig. 13 b). Fluxon modes realized at the end and the beginning of flux- flow branches for both sweeps 
are specified in Figs. 13 b) and c), respectively. Note that in both cases the maximum fluxon velocity is 
close to the lowest characteristic velocity ~ ci, i.e., splitting of c„ is not the cause of sub-branching of 
IVC's. Fig. 13 c) demonstrates that critical currents and return currents are also different for the two 
sweeps, i.e., Ic is multiple valued. As discussed in sec. 4.2. different /c's correspond to different modes, 
indicated in Fig. 13 c). 

Fig. 13 reproduces main features of experimental IVC's from Fig. 12. In Ref. it was shown 
that tiny sub-branches due to fluxon-antifluxon modes exist also at zero field. Switching between fluxon 
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modes/submodes and correlated fluctuations of flux-flow voltage would cause extremely wide radiation 
line width, consistent with numerical simulations [ p7[ and experimental radiation detection from HTSC 
SJJ's II . 

5 CONCLUSIONS 

In conclusion, single and multi-fluxon properties of layered superconductors and stacked Josephson 
junctions were analyzed. We have seen that behavior of SJJ's can be qualitatively different from that 
of single Josephson junctions or type-II superconductors. To a large extent, anomalous properties are 
due to unusual fluxon shape and existence of multiple quasi-equilibrium fluxon modes in long, strongly 
coupled SJJ's. 

The single fluxon in the stack of N junctions is described by N characteristic lengths. With increasing 
velocity, magnetic induction in a fluxon can change sign. For large N there is no Lorentz contraction of 
the fluxon at it — > ci . This results in absence of velocity matching behavior and a possibility of " super- 
relativistic" fluxon motion with u > ci. Such motion is accompanied by generation of Josephson plasma 
waves moving along with the fluxon (" Cherenkov" radiation) due to degeneration of fluxon components 
with c„i < u. 

Due to existence of multiple quasi-equilibrium fluxon modes/submodes the state of the stack is not 
well defined for given external conditions. It can only be described statistically with a certain probability 
of being in any of quasi-equilibrium states. This causes complicated and anomalous behavior of SJJ's: (i) 
Characteristics of SJJ's exhibit strong fluctuations, (ii) The c-axis critical current is multiple valued. The 
probability distribution P{Ic) has multiple maxima and Ic{T,H) consist of multiple sub-branches, (iii) 
The Ic{H) patterns are aperiodic, (iv) Flux-flow IVC's consist of multiple closely spaced sub-branches. 
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